set.seed(2423045)
rm(list=ls())
library(Matching)

source("AutoCluster3.R")

load("rep.sourcecode.re74psid1.RData")

X <-   as.matrix(cbind( foo$age,
                        foo$educ,
                        foo$black,
                        foo$hispan,
                        foo$married,
                        foo$nodegree,
                        I(foo$re74/1000),
                        #I(as.real(foo$re74 == 0)),
                        I(foo$re75/1000)))# ,
                        #I(as.real(foo$re75 == 0))))
                 
BalanceMat <- as.matrix(cbind(
                     foo$age,
                     foo$educ,
                     foo$black,
                     foo$hispan,
                     foo$married,
                     foo$nodegree,
                     I(foo$re74/1000),
                     I(foo$re75/1000),
                     I((foo$re74/1000)^2),
                     I((foo$re75/1000)^2),
                         
                     I((foo$age)^2),
                     I((foo$educ)^2),
                     I((foo$black)^2),
                     I((foo$hispan)^2),
                     I((foo$married)^2),
                     I((foo$nodegree)^2),
                         
                     I(foo$age*foo$educ),
                     I(foo$age*foo$black),
                     I(foo$age*foo$hispan),
                     I(foo$age*foo$married),    
                     I(foo$age*foo$nodegree),
                     I(foo$age*(foo$re74/1000)),
                     I(foo$age*(foo$re75/1000)),

                     I(foo$educ*foo$black),
                     I(foo$educ*foo$hispan),
                     I(foo$educ*foo$married),    
                     I(foo$educ*foo$nodegree),
                     I(foo$educ*(foo$re74/1000)),
                     I(foo$educ*(foo$re75/1000)),

#                     I(foo$black*foo$hispan), This is always zero!
                     I(foo$black*foo$married),
                     I(foo$black*foo$nodegree),
                     I(foo$black*(foo$re74/1000)),
                     I(foo$black*(foo$re75/1000)),

                     I(foo$hispan*foo$married),
                     I(foo$hispan*foo$nodegree),
                     I(foo$hispan*(foo$re74/1000)),
                     I(foo$hispan*(foo$re75/1000)),

                     I(foo$married*foo$nodegree),
                     I(foo$married*(foo$re74/1000)),
                     I(foo$married*(foo$re75/1000)),
                                                
                     I(foo$nodegree*(foo$re74/1000)),
                     I(foo$nodegree*(foo$re75/1000)),                        

                     I((foo$re74/1000)*(foo$re75/1000))                                                 
                         ) )


# From Review of Econ and Statistics
model.A = foo$treat~I(foo$age) + I(foo$age^2) + I(foo$age^3) + I(foo$educ) + I(foo$educ^2) + I(foo$married) + I(foo$nodegree) + I(foo$black) + I(foo$hispan) + I(foo$re74) + I(foo$re75) + I(foo$re74 == 0) + I(foo$re75 == 0) + I(I(foo$re74)*I(foo$educ))

pscores.A <- glm(model.A, family = binomial)

X2 <- cbind(X,BalanceMat[,9:ncol(BalanceMat)]) #first 8 spots of BalanceMat are already in X
for (i in 1:ncol(X2))
  {
    lm1 = lm(X2[,i]~pscores.A$linear.pred)
    X2[,i] = lm1$residual
   }

# VARIABLES WE MATCH ON BEGIN WITH PROPENSITY SCORE
orthoX2.plus.pscore = cbind(pscores.A$linear.pred, X2)
orthoX2.plus.pscore[,1] <- orthoX2.plus.pscore[,1] -  mean(orthoX2.plus.pscore[,1])


# NORMALIZE ALL X COVARS BY STANDARD DEVIATION
#for (i in 1:(ncol(X)+1))
for (i in 1:ncol(X2))
  {
    orthoX2.plus.pscore[,i] <-
      orthoX2.plus.pscore[,i]/sqrt(var(orthoX2.plus.pscore[,i]))
  }

#sv from psid1.re74.gm1.triton1.Rout
sv <- c(2.955110e+02, 1.426625e+01, 8.220941e-03, 8.967065e+02, 1.502436e+01,
        2.287152e-01, 2.816118e+02, 1.344551e+01, 5.291986e+02)
sv <- c(sv,rep(0,ncol(orthoX2.plus.pscore)-length(sv)))

#from psid1.re74.gm4.Rout, GENERATION: 29
sv <- c(2.955110e+02, 1.426625e+01, 8.220941e-03, 8.967065e+02,
        1.502436e+01, 2.287152e-01, 5.079356e+02, 1.344551e+01, 5.291986e+02,
        3.594517e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00,
        0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00,
        0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00,
        0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00,
        0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00,
        4.841528e+02, 0.000000e+00, 4.084006e+01, 0.000000e+00, 1.781796e+00,
        0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00)

#orthoX2.plus.pscore <- orthoX2.plus.pscore[,1:(8+cnew)]

for (s in 1:100)
  {
    cat("*S:",s,"**************************************************************\n")
    cl <- NCPUS()
    GM.out <-  GenMatch(Tr = foo$treat,
                        X = orthoX2.plus.pscore,
                        starting.values=sv,
                        BalanceMatrix = BalanceMat,
                        pop.size = 10000,
                        max.generations = 500,
	                hard.generation.limit=TRUE,
                        wait.generations = 50,
                        cluster=cl)
    stopCluster(cl)

    sv=diag(GM.out$Weight.matrix)
    cat("Starting Values:\n")
    write.csv(sv,eol=", ", row.names=FALSE)

    mout <- Match(Y=foo$re78,
                  Tr = foo$treat,
                  X = orthoX2.plus.pscore,
                  Weight.matrix=GM.out)
    summary(mout)
  }


